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1 Introduction 



High energy electron-positron or Bhabha scattering pQ is among the most important and care- 
fully studied processes in particle physics. It provides a very efficient tool for luminosity de- 
termination at electron-positron colliders and thus mediates the process of extracting physical 
information from the raw experimental data [2J. The small-angle Bhabha scattering is par- 
ticularly effective as a luminosity monitor at high-energy colliders]^ The large-angle Bhabha 
scattering is used to measure the luminosity at colliders operating at the center-of-mass energy, 
y/s, of a few GeV, such as BABAR/PEP-II, BELLE/KEKB, BES/BEPC, KLOE/DA$NE, 
and CMD, SND/VEPP-2M [5]Jf| Moreover, it will be also used to disentangle the luminosity 
spectrum at the ILC [6fT] . Bhabha scattering involves stable charged leptons both in the initial 
and the final states and, therefore, it can be measured experimentally with very high preci- 
sion. At LEP, the experimental error in the luminosity measurement has been reduced to 0.4 
permille [8] and it is expected to be even smaller at the ILC: the goal of the TESLA forward 
calorimeter collaboration is to reach the experimental accuracy of 0.1 permille in the first year 
of run [9]. Finally, at the low-energy accelerators DA$NE and VEPP-2M the cross section of 
the large-angle scattering is measured with the accuracy of about 1 permille pTOlTTT] . In the 
phenomenologically most interesting cases of low energy or small angle scattering, the Bhabha 
cross section is QED dominated, with the electroweak and hadronic effects being strongly sup- 
pressed. Therefore, it can be reliably computed in perturbative QED, with the accuracy limited 
only by uncalculated high order corrections. 

These properties make in such a way that Bhabha scattering is an ideal "standard candle" 
for electron-positron colliders. 

Realistic simulations of the Bhabha events, which take into account the detector geometry 
and experimental cuts, are performed by means of sophisticated Monte Carlo generators, such as 
BHLUMI P2], BABAYAGA [5lfT8] . BHAGENF [33], BHWIDE [T5], MCGPJ [TB], and SABSPV 
[TT] . To match the experimental needs, the two- loop QED corrections must be included into the 
theoretical analysis and incorporated into the event generators. Since the theoretical accuracy 
directly affects the luminosity determination and may jeopardize the high-precision physics 
program at electron-positron colliders, remarkable efforts were devoted to the study of the 
radiative corrections. The one-loop corrections have been known in the full electroweak theory 
for a long time [IE]. The two-loop electroweak corrections are still elusive. However, recently the 
calculation of the two-loop QED corrections was completed. These corrections can be divided 
into three main categories: (i) the pure photonic corrections, (ii) the corrections involving the 
electron vacuum polarization, i.e. with at least one closed electron loop, and (iii) the corrections 
involving the vacuum polarization by heavy flavors (leptons or quarks). The first results for 
the photonic corrections were obtained in the limit of small scattering angles pT9|2QII21|l22| . in 
the massless electron approximation [23J, and for the terms enhanced by powers of the large 
logarithm hi(s/m 2 e ) [23]. Finally, the photonic corrections to the differential cross section were 
obtained in [25J in the leading order of the small electron mass expansion through the infrared 
matching to the massless approximation. This result is sufficient for all phenomenological 
applications at present and future colliders [26] and was recently confirmed within a slightly 

4 At LEP, the luminometers were located at an angle between 1.4° and 2.9°. At the future International 
Linear Collider (ILC), they will be placed between 0.7° and 2.3° [3]. 

5 For example, in KLOE experiment the luminosity measurement is based on the events with scattering angles 
between 55° and 125° 0]. 
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different framework [27] (see also [28])|§ The corrections involving a closed electron loop were 
obtained in [35|36|l37j by direct diagrammatic calculation, retaining the full dependence on 
m e . The calculation was performed by using the Laporta algorithm [38] for the reduction of 
the Feynman diagrams to the master integrals (Mis) [39], which were subsequently evaluated 
[30|31f32|33|34] by means of the differential equation method [10]. The result was obtained 
in analytical form in terms of harmonic polylogarithms [Hll42|43!fl3] . The corrections due 
to the vacuum polarization by heavy fermion of mass > m e were first evaluated in the 
limit mj <C s,t,u by two different methods [2711^5] . The calculation of [27] is based on the 
expansion in the electron mass within the effective theory approach, while the calculation of 
[15] is diagrammatic and based on the reduction to the Mis evaluated in the asymptotic regime 
[46] . The condition m'j <C s,t,u, however, does not hold for r-lepton, c- and 6-quarks in the 
practically interesting energy range of about a few GeV, as well as for the top quark at typical 
ILC energies 500 GeV ^ y/s <^ 1000 GeV. In a recent letter jUJ we announced the result for 
the two-loop heavy-flavor contribution which is valid for any ratio of the heavy fermion mass 
to the Mandelstam invariants, provided s,t, u 3> m^. The calculation was performed in the 
small electron mass limit. We used the general theory of infrared and collinear divergencies 
to separate the singular dependence of the corrections on the vanishing electron mass. The 
most difficult part of the calculation was then carried out with a strictly massless electron. 
This critically reduced the complexity of the problem and made it solvable by the method of 
p|36f37] fl ' 

In this paper we provide a detailed account of our calculation [17] and we present the 
complete analytical result for the correction to the Bhabha cross section. The paper is organized 
as follows. In Section [2] we introduce our notations and conventions. In Section [3] we discuss the 
infrared and collinear structure of the corrections and outline the strategy of the calculation. 
In Section H] we describe the calculation of two previously unknown four-point two-loop master 
integrals. In Section we present the analytical result for the correction to the cross section. 
The numerical analysis is given in Section Section [7J contains our conclusions. Some technical 
aspects of the calculation including the auxiliary functions, generalized harmonic polylogarithms 
(GHPLs), and the asymptotic behavior of the corrections are discussed in the appendices. 

2 Notation and Conventions 

In this Section we briefly summarize our notation and conventions which follow [36 37] . We 
consider the photon mediated process 

e~(Pi) + e + (p 2 ) -> e~(p 3 ) + e + (p A ) , (1) 

where p\ = —m? e . In the following, we will neglect the electron mass, which is much smaller 
than any of the other mass scales involved in the problem. m e will be set consistently to zero 
everywhere, except where it acts as a regulator for the collinear singularities. The kinematics 

6 The full dependence of the pure photonic corrections on the electron mass m e is not known at the moment. 
The corresponding calculation involves the two-loop box diagrams with three scales: s, t and m e , which are not 
yet available, though many relevant results have been already obtained |29|30|31|32l33|34j . 

7 When this work was in preparation a numerical result for the two-loop heavy-flavor contribution was ob- 
tained by means of the dispersion relation approach [48j . 
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of the process is described in terms of the Mandelstam invariants s, t and u: 



s 
t 

u 



-P 2 = -(p 1+ p 2 ) 2 = 4E 2 , 

-Q 2 = -(pi-P3) 2 = -4^ 2 sin 2 ^ 

-R 2 = -( Pl -p 4 ) 2 = -AE 2 cos 2 ^ 



(2) 
(3) 

(4) 



where s + t + u = 0, E is the particle energy in the center-of-mass frame, and 9 is the scattering 
angle. The Bhabha scattering differential cross section is given by a series in the fine-structure 
constant a: 

da da$ /a\ do\ , / a \ 2 da^ 

dri ~ In + W ) In ' 



71 



dn 



+ O (a 3 ) 



(5) 



where 



da 

~dn 



a 
s 



st + — +f 



+ 



t 2 



t 2 , 
st + - + s' 



+ - s 

st 



+ 0(m 2 Js 



(6) 



is the Born cross section. In this paper, we consider only the radiative corrections dai/dVt 
that involve the vacuum polarization by heavy fermions. The first order correction to the 
cross section comes from the interference of the diagrams (a) and (b) in Fig. [U with the Born 
amplitude and it reads 



ail s J 



— [st + — + r 2 Rell 



r(u.o), 



s) + - [st + - + s 2 211 



r(li,0) 



+-(s + t) 2 (ReU^ \s) + U^°\t) 
st V 



+ 0{m 2 Js) . 



(7) 



The expression of the one- loop vacuum polarization functions IIq 1 '' ^ is given in Appendix [B} Qf 
is the electric charge of the heavy fermion, N c (number of colors) is equal to 1 for leptons and 
3 for quarks, and we adopt the on-shell scheme for the renormalization of a and of the fermion 
mass. Note that Eq. is infrared finite and has a regular behavior in the small electron mass 
limit. 



3 Structure of the Second-Order Corrections and Cal- 
culation Method 

The two-loop virtual corrections are infrared divergent. These soft divergencies are canceled 
in the inclusive cross section when one adds the photonic bremsstrahlung [19] . We regulate all 
the soft divergencies by dimensional regularization in D space-time dimensions. The standard 
approach to deal with the bremsstrahlung is to split it into a soft part, which accounts for 
the emission of the photons with the energy below a given cut-off ui m e , and a hard part 
corresponding to the emission of the photons with the energy larger than u. The infrared 
finite hard part is then computed numerically using Monte-Carlo methods with physical cuts 
dictated by the experimental setup. At the same time, the soft part is computed analytically 
and combined with the virtual corrections ensuring the cancellation of the soft divergencies in 
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(a) (b) 




(g) (h) (i) (j) 

Figure 1: One-loop diagrams. 

Eq. ([5]). Thus we consider the second-order contribution to the cross section given by the sum 
of two terms: 

dao daY da£ , , 

dn dn dn ' 1 ; 

which correspond to the two-loop virtual correction^] and the one-loop correction to the single 
soft photon emission which factorizes into the product of the first-order contributions [50J. 

The calculation of the virtual corrections is a highly nontrivial problem since in principle it 
involves the two-loop box diagrams depending on four mass scales: s, t, mf, and m e . These 
diagrams are beyond the reach of the available calculational techniques. However, in practice the 
electron mass is much smaller than any of the other mass scales involved in the problem and the 
calculation can be significantly simplified by exploiting this scale hierarchy. The small electron 
mass limit is not trivial because a finite electron mass regulates the collinear divergencies, 
giving a logarithmic dependence of the second order correction on m e . One way to perform a 
systematic expansion in the small electron mass is to use the expansion by regions approach 
[5Tj . However, if we are interested only in the leading order term in rn^/s, the problem can be 
solved in an elegant way without the expansion of the individual diagrams [37]. The main idea 
of the method is to use the general theory of collinear divergencies to identify a set of simple 
diagrams responsible for the singular behavior of the corrections on the electron mass. Then, 
we compute the remaining corrections with a strictly massless electron, effectively removing 
one mass scale from the most complicated part of the calculation. 

Let us describe the approach of [37] in more detail. The second order contribution to the 
cross section can be split in the sum of two terms according to the asymptotic dependence on 

8 We do not consider the trivial correction given by two heavy-fermion loop insertions which are are usually 
treated by Dyson resummation. 
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da 2 

In 



ln ( ) + C + frails) 



mi 



(2) 



dao 
dVL 



(9) 



The logarithmic term in Eq. is a remnant of the collinear divergence regulated by the 



electron mass. The quantities 5[ and 5 in Eq. ([9]) depend on s, t, and nif only. The 
collinear divergencies, and hence the singular dependence of the corrections on m e , have a 
peculiar structure which was extensively studied in the context of QCD. In particular, in a 
physical (Coulomb or axial) gauge the collinear divergencies factorize and can be absorbed in 
the external field renormalization [52]. Due to the factorization, the singular dependence on 
m e is the same for the Bhabha amplitude and the square of the vector form factor [25] . This 
attributes the total logarithmic corrections to the two-loop Bhabha scattering amplitude to 
the one-particle reducible diagrams (s)-(v) and the one-particle irreducible diagrams (g)-(j) 
of Fig. HI Moreover, due to the on-shell renormalization condition, the vacuum polarization 
does not change the photon propagator near the mass shell where the collinear divergencies 
are located. As a result, the irreducible diagrams are infrared finite even for m e = and the 
singular terms are entirely contained in the reducible diagrams. In calculating the cross section 
one has to take into account also the contributions coming from the interference of the one-loop 
corrections to the amplitude and the soft emission. Both contributions have a factorized form 
and can be easily evaluated for m e ^ 0. Thus, it is straightforward to obtain the coefficient of 
the logarithmic term in Eq. (Q, which reads 



(2) 



2 In 



4u x 



d<j\ 
dao 



(10) 



At the same time the sum of the remaining two-loop one-particle irreducible diagrams has a 
regular behavior in the small electron mass limit and can be computed with m e = 0. The 
two-loop vacuum polarization given by the diagrams (a)-(f) in Fig. H] does not develop collinear 
singularities, because the corresponding photon is far off-shell. Hence, the sum of the two-loop 
box diagrams (k)-(r) in Fig. H] is free of collinear divergencies as well. Let us emphasize that 
this property in general holds only for the sum of the diagrams. The individual diagrams 
computed in a covariant gauge do exhibit the collinear divergencies for m e = 0. This, however, 
does not pose any additional problem since we work in dimensional regularization. In this case, 
the collinear divergencies show up as extra poles in (D — 4), which are not related to the soft 
emission and disappear in the sum of the one-particle irreducible diagrams. The cancellation of 
the collinear singularities of the box diagrams in the Feynman gauge is schematically shown in 
Figs. [2] and [3] for the one- and two-loop cases, respectively. Let us demonstrate this cancellation 
explicitly in the case of the one-loop graphs. Each one-loop box diagram for m e = exhibits 
a double pole in (D — 4), arising from the overlapping of soft and collinear divergencies. In 
particular, for the diagram (g) in Fig. [1] one finds 



da{ 



v 



dn 



a 



1 



(g) 



s {D-Af 



Bp 



-2), 



S,t) 



0(l/(£>-4)) 



while for the diagram (i) one obtains 



da\ 



dn 



a 



(i) 



s (D-Af 



m 



f d(U-2), 



7Tb 

u,t)--+B^->(u.1) 



,(li,-2) , 



+ 0(1/(D - 4)) . (12) 
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= Free of collinear poles 



Figure 2: Cancellation of the collinear poles among one-loop box diagrams calculated by setting 
m e = from the start. 




= Free of collinear poles 



Figure 3: Cancellation of the collinear poles among two-loop box diagrams calculated by setting 
m e = from the start. 

The explicit expressions of the auxiliary functions B^ 11 ' ^ (z = 1,2,3) for x > are collected 
in Appendix [HI It is easy to check that 

B[ ll '- 2 \s,t)+Bi ll '- 2 \u,t) = 0, 

B^- 2 \s,t)-B^- 2) (u,t) = 0, (13) 

so that the double pole disappears in the sum of the diagrams. The residual single pole in (D— 4) 
is of soft nature and it is canceled after adding the soft-photon emission. The cancellation of 
the collinear poles of the two-loop box diagrams is completely analogous to the one-loop case. 

Since the sum of the box diagrams has a smooth limit m e — > 0, the result does not depend 
on whether this limit or the limit e — > is taken first. In other words, the absence of collinear 
divergencies makes in such a way that the expression of the sum of the box diagrams cannot 
depend on which collinear regulator (electron mass or dimensional regularization) is employed 
in the calculation. All the "true" two-loop diagrams contribute only to the non-logarithmic 
term in Eq. (Q and, thus, can be evaluated for m e = 0. The two- loop problem with massless 
electron falls in the same complexity class as the one considered in [35|36|37j and can be solved 
by similar approach. In the reduction of the two-loop box diagrams, however, two completely 
new Mis appear. The calculation of these Mis is described in the next section. 

4 The Master Integrals 

The two-loop heavy-fermion correction to the Bhabha scattering amplitude is given by the 
Feynman diagrams shown in Fig. HI We express the square modulus of the amplitude in terms 
of scalar integrals. The ultraviolet, soft, and collinear divergencies of the integrals are treated 
by dimensional regularization. By means of the Laporta algorithm [38] the scalar integrals are 
reduced to six Mis diagrammatically shown in Fig. [5j Four of them, Fig. [5] (c)-(f ), were already 
known [3Tf53] . The integrals Fig. [5] (a) and (b) represent the main computational result of the 
present paper. Below we describe their calculation. 



6 



(pa ■ fa) 





(a) (b) (c) 

(j>3 ■ k 2 ) 





(d) (e) (f) 

Figure 5: The set of 6 two-loop Master Integrals involved in the calculation. 

For the pair of Mis of the box topology we choose the integrals Fig. (a) and (b) with the 
following momentum routing: 



M^D^lP^Q 2 ) = () — ^— — , (14) 



M 2 (D,m},P 2 ,Q 2 ) = j Q (pa fc 2 ) = J ® D h® D k 2 - - , (15) 



V{D 3 V 4 V 5 V 6 

where Pi = fcf , £> 2 = (p a - £;i) 2 , P 3 = (?2 + ^i) 2 , P 4 = ^ + m 2 , P 5 = (pi -p 3 - h + /c 2 ) 2 + 
and where the integration measure is defined as 



m 



f 



C(D) is a function of the space-time dimension D: 



C(L>) = (4tt)— r(3--j , (17) 

with C(4) = 1. In Eq. ffl6|) /i stands for the 't Hooft scale of dimensional regularization and we 
set fi = rrif in the rest of the paper. The integration measure in Eq. ffTBT) is chosen in such a 
way that the one-loop massive tadpole becomes 

v D k T ^-^= ,„ m L .. . (18) 



P + m 2 p-2)(D-4)" 

The Mis M\ and M 2 satisfy two systems of first-order linear differential equations [101 in the 
external kinematical invariants s and t. With our choice of Mis, the s-system is completely 
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decoupled, i.e. every MI satisfies a single first-order linear differential equation. By contrast, 
the t-system consists of two coupled equations and, therefore, is equivalent to a second-order 
linear differential equation for one of the MlsE The solution of the system requires two initial 
conditions. Our Mis are functions of ml, s and t. They are regular for t — > and it is possible 
to determine the initial conditions at t = for any value of s. Therefore, it is more convenient 
to solve the system of the differential equations in t then the decoupled first-order equations in 
s. We define the following dimensionless variables: 



x 



p2 

2 

m 



In terms of these variables the i-system takes the form: 



2 ' 



y 



2 

m 



m. 



(19) 







~D-5 


dy 




[ 2y 


dM 2 






dy 


2(y + A) 



D-A 



2(y+A) 2(y+x 
D — A 



Mi - 



I 



f 

- + 

y y + x 



3D-10 
2m) 

3D- 



1 

y 

10 



M 2 + n 1 {D,m 2 f ,x,y), (20) 



2(y + 4)J 



y+A 

M 2 + n 2 (D,m 2 f ,x,y). 



(21) 



Note that the second MI with the scalar product on the numerator is dimensionless, while the 
first one has the mass dimension of mj 2 . The functions £li(D, m 2 , x, y) are linear combinations 
of the Mis (c)-(f ) of Fig. El the product of a massless one-loop bubble and a massive tadpole, 
and the product of two tadpoles: 



Q 1 



x(x — 4) 



y(A + y)(x + y) 
1 




2mj(D - A)y 2 (y + A) 2 (y + x) 



-560 + 384£> - 64L> 2 + 200x - lAODx + 2AD 2 x 




+72y - 76Dy + 16D 2 y + lODxy - 3D 2 xy + A8y 2 - 40/V + 8D 2 y 2 ) 
12(D - 2) (-28 + 8D + lOx - 3Dx - 6y + 2Dy) (p 3 . fc 2 ) 

mj(D - A)y 2 (y + A) 2 (y + x) 

2(D-2)(D-3) r , 
mj(D — A)xy(y + 4) _ V>" X ^ 

(D - 2) 




8m 6 f (D -5)(D- A){D - 3)y 2 (y + A) 2 (y + x 



■(5600 - A960D + 1A08D 2 - \28D C 



-2000s + \8mDx - 520D 2 x + A8D 3 x + Utey - 1096Dy + 328D 2 y - 32D 3 y 



-296xy + 30ADxy - 98D 2 xy + IQD^xy - 2Ay 2 + 2$Dy 2 



AD 2 y 2 



2Axy z 



Q 9 



+26Dxy 2 

(D-A)x 
8y(y + x) 



9D 2 xy 2 + D 3 xy 2 ) QQ 

m 2 x(x — A) 
A(y + A)(y + x) 



(22) 




+ 




-360x + 22ADx - 32D 2 x - 480y + 3l2Dy - A8D 2 y 



Amj(D - A)y 2 (y + A) 2 (y + x) 

9 As it is shown below, the second-order equation in t turns out to be particularly simple because the system 
of the first-order equations in s is decoupled. 
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+16xy - 3ADxy + 8D 2 xy - 8y 2 - 20Dy 2 
+4% 3 - AADy 3 + 9D 2 y 3 ) 



QD 2 y 2 + 2Axy z - 21Dxy z + AD' z xy 2 




Q(D - 2)(18x-ADx + 2Ay-6Dy+Axy-Dxy + Wy 2 - 3Dy 2 ) 
mj(D - A)y 2 (y + A) 2 (y + x) 
(D-3)(D-2) 




(p3 ■ k 2 ) 



2m 2 (D-A)x(y + A) 



o 



D — 2 



■(7200x - 5920Dx + 1536D 2 x 



32mj(D -5)(D- A)(D - 3)y 2 {y + A) 2 (y + x 
-128D 3 x + 9600y - 8160Dy + 2208D 2 y - 192D 3 y + 1360x?/ - 1032Dxy 
+232D 2 xy - 16D 3 xy + 3424y 2 - 29UDy 2 + 808D 2 y 2 - 72D 3 y 2 - llQxy 2 
+200Dxy 2 - 12D 2 xy 2 + 8D 3 xy 2 + 16y 3 - 8Dy 3 - 24xy 3 + 26Dxy 3 - 9D 2 xy 3 
+D 3 xy 3 ) OO • ( 23 ) 

All the possible singularities of the integrals Mi and M2 are those appearing in their coefficients 
in Eqs. ( I20f21j) . Thus, the singularities are potentially located at y = 0, y = — 4 and y = —x. 
The point y = —4 is a singular point for the integrals. It corresponds to the three (two massive 
and one massless) particle cut in the t channel. By contrast, y = is a regular point. We can use 
this information in order to determine the initial conditions. In fact, multiplying Eqs. ( |20l I2TI) 
by y and taking into account that y dMi^/ dy\ y ^$ — > 0, in the limit ?/->0we find 



M 1 (D,x,y = 0) 



G{0,x) 



+ 



C(2)-G(0,0,x) 



lm 2 x 



(D — 4) 



2m 2 f x (D - 4) 3 4m 2 f x (D - 4)< 

J - 2C(3) + (C(2) - 4)G(0, x) - G(0, 0, 0, x) + G(jjl, (i, 0, x)] 



16m 2 x 



A- x 



8m 2 f x\J 'x(4 



-.G(n,0,x)+O(D-4) 



(24) 



x 



1 



M 2 (D } x } y = 0) = — 



1 



16 (D — A) 2 



(2-G(0,x)) 
32 



1 



1 



\llx-({2)x-5xG{0,x) 



(D - 4) 64x ' 

+xG(0, 0, x)+2G(n, /i, 0, x)}+ ~ X G(n, 0, x) + 0(D - 4) , (25) 

64^x{4— x) 

where the functions G are generalized harmonic polylogarithms (GHPLs) [SI] described in 
Appendix E The system of Eqs. (120112 ip is equivalent to a second-order linear differential 
equation for one of the Mis. For the integral Mi, after trivial manipulations, we find 



d 2 M 1 
~dy~ 2 



+ 



3 

D - 6 



3D - 13 D - 4 



2(2/ + 4) 



y + x 
3D - 



dMi 
dy 



-{D-4) 



D-5 



13 



A(y + x) 2 
+tt(D,m 2 f ,x,y) 



4(x — 4) J y + x 



4y2 

3D - 17 



16 



+ 



3D - 17 
16 

3D - 13 



Ax J y 



A(x-A) J y + 4 



Mi 



0. 



(26) 



The function Q(D,m 2 ,x,y) contains the Mis (c)-(f) of Fig. some products of one-loop 
integrals and the product of two tadpoles. 
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Let us discuss briefly the structure of the Eq. (1261) . Note that we are not interested in a 
solution of the differential equation valid for arbitrary value of the space-time parameter D. 
Instead, we look for a solution in the form of a Laurent series in (D — 4). Since the coefficient 
of Mi is proportional to (D — 4), Eq. (1261) at each order in (D — 4) is reduced to a first-order 
linear differential equation for the derivative M[ = dM\jdy: 



N 



M[(D,m},x,y) = ^ M' li (D — 4)* + 0(D — A) N+1 



i=-3 



dM[ t 
dy 



3 1 

+ 



2y 2(y + A) 



(27) 
(28) 



The functions ^/ i (m'j,x,y) in Eq. (T2"gj) are defined as follows {M 1 j = for j < —3): 



-Qi(D,m 2 f ,x,y) + 



2(y + A) y + x. 



Ki-i 



i 



1 



3 

— + 



3 5 
4x 16 J y 
1 



+ 
+ 



2(y + x) 2 \Ax A{x-A)J y + x \ 16 4(x - 4) / y + 4 
"13 1 3 1 

Ay 2 16y A(y + x) 2 A (x - 4) y + x 



+ 



16 A(x -A) J y + A 



1 



i-2 



(29) 



Note that ^i(m 2 ,x,y) contains M[^ i _ 1 , Mi^-i and Mi i j_ 2 - Moreover, fj(mj,3;,t/) contains the 
i-th order term of the Laurent expansion of Q(D, m 2 , x, y): 



o 



tt(D,m 2 f ,x,y) 



Vi(D, m}, x, y) (D - 4) i + 0(D - 4) 

i=-2 

2a; 2 + 10a;?/ - x 2 y + Ay 2 1 
Am 2 f xy 2 (A + y)(x + y) 2 (D - A) 2 



8m 2 y 2 (x + y) 2 
2x + 6y — 



3x + 20y — 2xy + y 2 2x 2 + 10a;?/ — x 2 y + 4y 5 



G(-fx,y) 



(4 + 2/) 
1 

. P^4) 



x(4 + y) 



G(0,x) 



16a; 2 + 40a;y - 4x 2 y + 2xy 2 - 2x 2 ((2) 



i r 

lQm 2 y 2 (x + y) 2 x(4 + y) \ 
-10xyC(2) + x 2 yC(2) - 4y 2 C(2) - a;(4a; + 16y - 3xy)G(0, x) + (2a; 2 
+10xy - x 2 y + 4y 2 )G(0, 0, x) + 2x(2x - 2y + xy)G(-(i, -//, y) 



xy(x-4) 2 x(y + 4) 

=G(yx, (J, x) 



(12a: + 24y-a;y + 2y 2 )G'(- Ai ,y) 



V^(4 - x) ' ' vW+4j 
-(2a; + 6y - xy)G(0, -/i, y) - 3(2a; + 6y - xy)G{-A, - M , y) 
+£>(£> -4). 



(30) 
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The formal solution of Eq. ( 1281) reads 

ML ' 



yVv(v + 4 ) 



/y 
r^Jr(r + 4)^ dr + J { 



(31) 



where Jj are integration constants, which are fixed by the regularity condition for the derivative 
M[ i at y = 0. Actually, in every order in — 4) we find Jj = 0. This means that we have to 
discard the homogeneous solution and keep only the particular solution of the inhomogeneous 
Eq. fl2HP. Then, we integrate Eq. (jglj) : 



Mi. 



M[ i (m 2 f , x, r) dr + . 



(32) 



Ji t i are again integration constants which are determined by the initial conditions Eq. (1241) . 
Once we have the master integral Mi, the calculation of M2 is straightforward. From Eq. (12 0j) . 
we express M2 in terms of Mi and dM\jdy. The final analytical expression for Mi and M2 
reads 



Mi(D, m% x;y)=J2 M^m}, x; y) (D - 4)* + 0(D - 4) 



i=-3 



where 

M l5 _ 3 

M l5 _ 2 
Mi,_i 



2m^x ' 



imjx 



2-G(0;x)- l±JL= G(-n;y) 



y(y + 4) 

4 + C(2) + 2G(0; a;) - G(0, 0; x) + 2G(-//, -fi; y) 



(33) 

(34) 
(35) 



y + 4 



M10 



(x-4)(y + 4) 



2G(-//; y) - 3G(-4, y) - G(0; x)G(- M ; y) 
G(/i, 0; x)G(-x - n; y) 



(36) 



16m^ a/x(4 — x)^Jy(y + 4) 
4 + y 



=[-G(-x,0,-//;y) + 6G(-/i, -//;y)-(4-C(2))G(-//;y) 

16mjxy / y(y + 4) 

+G(0; x) (-3G(-4, -//; y) + G(-x, -//; y) + 2G(-/i; y) - G(0, -//; y)) 
+6G(-4, -//; y) - 9G(-4, -4, -//; y) + G(0, 0, y) - G(0, 0; x)G(-//; y)] 

+ ^^[8 - 2C(2) - 2C(3) - 4G(-/i, -//; y) + 6G(-/i, -4, -//; y) + G(/i, /1, 0; x) 

-4G(0, -/i, -//; y) - (4 - C(2))G(0; x) + 2G(0; x)G(-/i, -//; y) + 2G(0, 0; x) 
-G(0,0,0;x)], (37) 



and 



M 2 (D, m 2 f , x;y)=J2 M ^A m % v) ( D ~ 4 )* + °( D ~ 4 ) 



(38) 



i=-3 
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where 

M 2j _ 3 

M 2 ,-2 

M2.-1 



M- 



2.0 



y_ 

8x 



16x 
1 

32x 
4 



x + 2/(4 - G(0; x)) - y) 



(39) 
(40) 



V^+4) 

C(2)y-5x-14y + 2(2+y+-)G(-/x, -/x; y) + (x+4y)G(0; x) - yG(0, 0; x) 

[(x + 4y)G(-/x; y) - 3yG(-4, -/x; y) - yG(0; x)G(-/x; y)] , (41) 

x - 4 



44-y 



32a; v ^T4) 

y(y + 4)(x-4) 



64 v /x(4-x) v /y(y + 4) 

y(y + 4) 



G(/x, 0; x)G(— x — /x; y) — 



64- v /x(4 - 

. G f(_ a . ,- A i;y) + (C(2)-14-— )G(-Ai;y) 

64a; a/?/ (y 4-4) L V y/ 

+ (l2 + — )g(-4, -/x; y) + 6G(-/x, -/x, -/x; y) - 9G(-4, -4, -/x; y) 

V y / 



G(mu } 0; x) 



x 



+-G(0, -/x; y) + G(0, 0, -//; y) + G(0; x)(G(-x, -//; y) - 3G(-4, -/x; y) 

y 

+4G(-/i; y) - G(0, -/x; y)) - G(0, 0; x)G(-/x; y) 1 ' 



64a; 



19x + 46y-C(2)x 



-4C(2)y - 2C(3)y - 2 (e + 3^ + x + 4y) G(-/x, — /x; y) 

+6 2 + - + y G(-/i, -4, -/x; y) + — G(-/x, 0, -/x; y) + (2 + y)G(/x, /x, 0; x) 
v 1/ / y 

-2 (4 + - + 2y)G(0, -/x, -//; y) + 2(2 + y)G(0; x)G(-/x, -/x; y) 
+(C(2)y - 5x - 14y)G(0; x) + (x + 4y)G(0, 0; x) - yG(0, 0, 0; x) 



(42) 



The result is expressed in terms of GHPLs of two variables, x and y. More details on these 
functions can be found in Appendix [A] We checked that Eqs. (I33|I38 j) do satisfy the system of 
linear differential equations in s, and the initial condition of Eq. (I25I) for M2 is recovered in the 
limit x — > 0. 



5 The Two-Loop Heavy-Flavor Correction 



In this Section we present the analytical result for the contribution of the different classes of 
the two- loop diagrams shown in Fig. H] and the corresponding soft-photon emission contribution 
to the differential cross section. We keep the notations as close as possible to [361137] . The 
ultraviolet divergencies are renormalized in the on-shell scheme. 

It is convenient to split the two-loop virtual correction in Eq. (jSJ) into the sum of five terms 



da\ 



da\ 



dVL dVL 



+ 



da\ 



(21, s) dVt 



+ 



daX 



(2i,v) dQ 



+ 



da\ 



(21, B) dVt 



+ 



(2l,R) dVt 



+ 



da\ 



(s,v) dQ 



(S,B) 



(43) 
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(a) (b) (c) (d) 




(e) (f) (g) (h) 

Figure 6: Diagrams contributing to the real corrections to the NNLO heavy flavor cross section. 

which correspond to the contribution of the two-loop self-energy diagrams, two-loop vertex 
diagrams, two-loop box diagrams, two-loop reducible diagrams, and to the interference of one- 
loop vertex and one-loop box diagrams with the one-loop self-energy diagrams, respectively. 
We drop the list of arguments of the various contributions to the cross section. All the terms 
in Eq. (14*3]) depend on s, t, and rrif. As it is explained in Section [3], only the fourth and fifth 
terms on the r. h. s. of Eq. (j43p depend logarithmically on the electron mass m e . 



5.1 Two-Loop Vacuum Polarization Corrections 

The contribution of the diagrams Fig. H] (a)-(f) can be obtained by replacing ITq 1 ^ with the 
two- loop vacuum polarization function 1I5 in Eq. (171): 



a 



(2l,S) 



QjN c 



st + — +f 



2ReII 



(2/), 



S) + 



t 2 



st + - + s z 



2Uf\t) 



+ 



l(, + t) 2 (Renf) ( s )+nf ) (t) 



(44) 



The explicit expression of the renormalized function ITq 2 ^ in terms of GHPLs is given in Eq. (11621) 
of Appendix [Bj Note that the two-loop vacuum polarization corrections are proportional to 
the fourth power of the heavy fermion charge, while all the other corrections are proportional 
to Q). 



5.2 Two-Loop Vertex Corrections 

The two-loop vertex diagram are shown in Fig. H] (g)-(j). These diagrams are infrared finite 
and can be evaluated for m e = 0. The analytical result for the two-loop vertex correction reads 



(2l,V) 



2—Q}N c 
s 3 



ivf >(m) + ivf >m) + 1 (vr\s,t )+ vr\t,s)) 



(45) 
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The functions v} 21 ^ and V r 2 vzt; in this equation are related to the two-loop Dirac form facto: 
Fi\—p 2 ) as follows: 

V} 2l \s,t) = Cljl (s,t)ReF^\t), 
V 2 (2l \s,t) = c 2A (s,t)ReFi 2l \t) , 

where the coefficients c^i(s,t) read 



'(») 



10 



(46) 
(47) 



ci,i {s, t) = (s + t) 2 , c 2 ,i (s, t) = 2 1st + s 2 + 



(48) 



The explicit expressions of the renormalized form factor is given in Eq. (11791) of Appendix [Bl 
5.3 Two-Loop Box Corrections 

There are eight two- loop box diagrams shown in Fig. [4] (k)-(r). If we transform the external 
momenta p^ «-> —p\ and P2 <-> — £>3 in the diagram (1), it becomes identical to the diagram (k). 
Since this transformation does not change the Mandelstam variables s and t, the contributions 
of diagrams (k) and (1) to the differential cross section are equal. The same is true for the pairs 
of diagrams: (m)-(n), (o)-(p), and (q)-(r). Thus the contribution of the two- loop box diagrams 
can be written as follows: 



da\ 



(21, B) 



2 2 

2—Q 2 f N c \^(ReB (2l \s,t) + ReBi 2l \t,s) + B\ 2l \u,t)-ReB {2l \u,s' 
4s 1 L s V 



m. 



+ -L[ReB {2l \s,t) + ReB\ ZL, (t,s) - B? l) (u, t) + ReSf >{u, s) 



>(20, 



>(20/ 



>(«)/ 



(49) 



where the overall minus sign is due to the closed fermion loop and the overall factor 2 reflects the 
identity of the diagrams discussed above. The Laurent expansion of the renormalized functions 
B^ reads 



Bf l \s,t) 



(D-A) 2 



B 



(21,-2) 



(s,t) + 



(D-4) 



B 



(21,-1) 



(s,t) + B^°\s,t) + 0((D-4)) . (50) 



The expressions of the coefficients Bf 1 '^ for j = —2,-1,0, are given in Eqs. (jl81H189"l) of 
Appendix [Bl As it was pointed out in Section [31 the double poles of the auxiliary functions 
Bi cancel in the sum of the planar and crossed box contributions. 



5.4 Two-Loop Reducible Corrections 

The two-loop reducible diagrams are shown in Fig. [4] (s)-(v). They reduce to the product of 
the one-loop vacuum polarization function and one-loop Dirac form factor and can be easily 
obtained from Eq. (43) of [36]. The result is of the following form: 



dn 



(2l,R) 



2—Q 2 f N c 
s J 



- 2 V 2 R (t, s) + ±V 2 R (s, t) + i (V*(s, t) + Vfit, a)) 
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The Pauli form factor vanishes in the limit m e 



t 2 
0. 



(51) 
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where the functions V i (s, i) have the Laurent expansion 



L-V^-^t) + Vf R '°\s,t) + o((D- 4)) 



(D-A) 



with 



V,; 



c J , 1 ( S ,t)Re(F 1 (1W) (t) nf' 0) (i) 
8 ,t) = cu(s,t) Re( F™(t) U^°\t) + if ^(i) U^\t) 
+c ij2 ( S ,i)Re(if Z '- 1) (i) II<"' 0) (t)). 

The coefficients ci.i(s, t) and C2,i(s, i) are given in Eq. ( 1481) . Moreover, we have: 

1, , , 1 



ci, 2 (s,t) = Ust + s 2 + t 2 ) 



c 2 , 2 (s,t) = ^t 2 . 



(52) 



(53) 



(54) 



(55) 



In Eqs. ( 15311541) the quantities fl^ (i = 0, 1) and F x ( 



(U,i) 



—1, 0) are the coefficients of 



the Laurent expansion of the one-loop vacuum polarization function and one-loop Dirac form 
factor given in Eqs. (11611 1 1 62[) and Eqs. (11641 II 65[) of Appendix [Bj respectively. Note that the 
form factor should be calculated by keeping a nonzero electron mass as collinear regulator. 



5.5 Interference of Two One-Loop Graphs 

Finally, we discuss the interference between the one-loop vacuum polarization diagrams of 
Fig. [1] (a) and (b), and the one-loop vertex and box diagrams of Fig. [1] (c)-(j). The calculation 
is straightforward because the one-loop vacuum polarization factorizes with respect to the tree- 
level amplitude and the interference term can be obtained from the one-loop vertex and box 
corrections to the cross section. For the vertex diagrams we obtain 



o 



(S,V) 



2—Q 2 f N c Re 



i^"'(M)(n«" , ( s ))" + i^"'( s , ( )(nJ">( ( ) 



+ iK< 1 "( s ,t)(n< 1 "w)* + ^r ( M)(<' , W 



where the functions (s,i) are defined as 

Vf>(M) = 
V™{s,t) = 
Their Laurent expansions read 



2 ( st + ^s 2 + l -A + l(D - 4) (si + s 2 + t 2 ) 



ReFl ll \t) 



■2[st + s 2 + U 2 \ + l -(D-A)t 2 



ReF} ll \t) 



i 



(D-A) 



t) + v} ll >°\s,t) + o((D-4)), 



(56) 



(57) 
(58) 



(59) 



16 



where 



(s,t) 
s,t) 



(60) 
(61) 



and with the coefficients Cij(s,t) given in Eqs. (I48|55l) . Thus, if for instance we consider the 
the first term of Eq. (1561) . we obtain 



+V^~ 1 \t,s){u^ 1 \s)y + 0[{D- 4)) . 



(62) 



Similar expressions hold for the other terms. 
For the box diagrams we obtain 



da\ 



dtt 



(S,B) As^ f 



Q)N c Re 



f - (b^ (s, t)+Bi ll) (t, s) + B£ l) (u, t)-Bi ll) (it, s)) (n 
+^(B^(s lt )+B[^s)-B^ 



(63) 



where the functions B^ (s,t) (i = 1,2,3) have the Laurent expansion: 



(1Z.-2) 



M) + 



-B 



s,t) + B^ \s,t) + O((D-4)) , (64) 



(D-A) 2 1 v ' y CD -4)' 

The explicit expressions of the coefficients of the Laurent expansion can be found in Appendix IB1 
Thus, for instance for the first term in Eq. f l63l) we obtain 



B?\s,t)(ur\s)y = ^^^(^^(nr^jV^rHM)^ 

+B^- 1 \s,t){u^ 1 \s)y + o((Z>-4)) , 
and similar expressions for the other terms. 



(iJ,o), 



(65) 



5.6 Soft Photon Emission 

Let us now discuss the calculation of the second term in Eq. (jSJ). We use the procedure applied 
in [37] to the case of the electron vacuum polarization. It is convenient to introduce the quantity 

da{°~ 4) 



dQ 



d<7 1 , 
(ii,s) d£l dVt 



where the first term on the right hand side (r.h.s.) is defined in Eq. flTJ) and the second term 
for m P = reads 



da[ 



(D-4) 



dn 



a 

(ll,S) s 



st + — + r 



2Reni 1 '' 1) (s) + 



t 2 



t 5 
St + — + s z 



2IT^ 1) 



(t) 



+±(. + tf (Renr\s) + nr\t)) + ReIV (s) 2 +n ^ (t) 



(67) 
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The contribution of the soft-photon emission is then given by 



7 7 = 1 



h m e ,m f ) 



Here, the infrared divergent quantities J\j (j — 1, • - - ,4) are defined as follows: 

Jy(pi ■ Pj, ml, m)) = €j (p! ■ pj) /^(p! • Pj, TTtg, mj) , 

where e.,- = +1 for j = 1, 4, = —1 for j = 2, 3, and 



hj(pi-Pj,m e ,m f ) 



in 



r (3 - f ) -K^-i)/2 A7X 2 J ko (p x • k)(pj ■ k) 



d u k 



1 



(6£ 



(69) 



(70) 



The calculation of the integrals Jy follows the procedure outlined in [55] and is described in 
detail in Appendix A of [37J. We need only the leading terms of the small electron mass 
expansion of these integrals, which to O ((D — 4)°) read 



hi 



'12 



2m 2 . 
In 



D-A 



mi 



m: 



+ In 



D -4 



+ In 



+ In 

Au 2 



-In 



ni 



f, 
+ In 



mi 



O(m e ) 



(71) 



In 2 



mi 



mi 



~ 2C(2) 



(72) 



'13 



In - 



-2C(2) 



mj 



O 



D-A 
,:2 



+ In 



Au z 



+ In 



In 2 



ml 



mi 



- Li 2 1 



(73) 



'ii 



In I - A 

mf 



-2C(2) 



+ '-z 



D-A 



mi 



In 



4u; s 



s 

9 

ml 



In I -3 



mj 



-In 2 — -Li 2 1 + - 



(74) 



Note that the term proportional to (D — 4) in Eq. (I6"6"|) gives a finite contribution to Eq. 
since contain an infrared pole. 



6 Numerical Analysis 

In this Section we consider the application of our result to the phenomenologically interesting 
cases relevant for physics at DA$NE and the ILC We provide a detailed account of the pure 
QED contribution, extend the analysis to the mixed QED-QCD corrections, and give an esti- 
mate of the hadronic vacuum polarization effect. All the terms involving the logarithm of the 
IR cut-off u, ln(4a; 2 /s), are excluded from the numerical estimates because the corresponding 
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u flO" 4 ) 


c (IO -4 ) 


r flO" 4 ) 


b (IO -4 ) 


ou 


1 1 .041UU4 


1 7079877 
L.I 3 1 ZO ( ( 


n nfi99fi77 


U.UZ04U1O 


U.UUlUozo 


ou 


lo.4U ( ooD 


9 99fi7Rt;/l 


n nsRi S7R 

U.UoOlo ( 


U.UoO ( UOo 


n nm a 1 s/i 

U.UU141o4 


7n° 

( u 


iy .400 i io 


9 R^nzLQ^n 
z.uou4you 


n i nsfii 9R 

U. lUoUlZU 


u.u4uoozy 




80° 


20.465455 


3.0655973 


0.1253094 


0.0540991 


0.0022442 


90° 


21.463240 


3.4581845 


0.1321857 


0.0576348 


0.0024428 


100° 


22.366427 


3.8070041 


0.1268594 


0.0560581 


0.0024304 


110° 


23.099679 


4.0922189 


0.1098317 


0.0495028 


0.0022024 


120° 


23.605216 


4.3030725 


0.0843311 


0.0392810 


0.0018086 


130° 


23.847394 


4.4392717 


0.0549436 


0.0273145 


0.0013297 



Table 1: The second-order electron, muon, c-quark, r-lepton, and b-quark QED contributions 
to the differential cross section of Bhabha scattering at ^/s = 1 GeV in units of 10~ 4 of the 
Born cross section. 



contribution critically depends on the event selection algorithm and cannot be unambiguously 
estimated without imposing specific cuts on the photon bremsstrahlung. The actual impact of 
the two-loop virtual corrections on the theoretical predictions can be determined only after the 
result of the paper is consistently implemented into the Monte-Carlo event generators. Never- 
theless, the above naive procedure can be used to get a rough estimate of the magnitude and 
the structure of the corrections. 

As a first application, we consider the Bhabha scattering at y 7 ? = 1 GeV. The latter is 
the value of the center-of-mass energy of the KLOE experiment at DA$NE, which plays a 
crucial role in the determination of the hadronic vacuum polarization contribution to the muon 
anomalous magnetic moment [TUfTI]. The anatomy of the heavy-flavor two- loop correction for 
this choice of the center-of-mass energy is shown in Fig. [?l where the case in which the heavy 
fermion is the r-lepton is considered. Each curve plotted in Fig. [7] represents a specific subset of 
the virtual corrections and the corresponding soft emission. The dominant contribution origi- 
nates from the two-loop reducible corrections and from the product of the one-loop corrections, 
which are considered together (see Sections 15.41 and 15. 5j) . This subset is numerically dominant 
because it contains all the large collinear logarithms ln(s/mg). The numerical values of the 
ratio of the second-order heavy-flavor corrections to the Born QED cross section (Eq. ([6])) are 
collected in Table [Tll"1 We separately consider the contributions of muon, r-lepton, c-quark 
and 6-quark. Since at KLOE one is particularly interested in the large angle scattering events, 
we considered the angular range between 50 and 130 degrees. For comparison we also give the 
value of the electron vacuum polarization contributiono The following input parameters are 
used: m e = 0.510998902 MeV, m M = 0.105658369 GeV, m T = 1.7 GeV, m c = 1.25 GeV, and 
nib = 4.7 GeV. The contributions of the r-lepton, c- and 6-quark are suppressed with respect to 
the muon at least by one order of magnitude. The total heavy-flavor contribution is dominated 
by the muon loop and it reaches 0.45 permille in magnitude at 6 ~ 140°. Note that the energy 

11 Note that the numerical evaluation of our analytic formulas is done with double Fortran precision but we 
do not present here all the available significant digits 

12 The electron contribution includes the logarithmic part of the soft-pair production. The logarithms of the 
soft-pair cut-off are excluded from the numerical estimates 26 . 
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Figure 7: Self-energy ("S"), vertex ("V"), reducible plus one-loop times one-loop ("R"), and 
box ("B") contributions to the two-loop r-lepton correction to the differential cross section of 
Bhabha scattering at yfs = 1 GeV. 



under consideration is sufficiently below the quarkonium threshold so that J Amj — s ^> Aq C d 

and the heavy quarks can be treated perturbatively. Moreover, for the heavy-quark vacuum 
polarization in two-loop approximation one has to take into account the first order corrections 
in the strong coupling constant a s due to a gluon exchange inside the quark loop. The resulting 
0(aa s ) correction to the Bhabha cross section can be obtained from the QED contribution, 
Eq. fj4"3]) . by adjusting the overall factor: 



do% Qcd c F a s (m 2 f ) do\ 
f 



dQ 



(2i,s) Q 2 f a dQ 



(75) 

(21, S) 



Cp = (N 2 — 1)/(2N C ) is the Casimir operator of the fundamental representation of the SU(N C ) 
color group, and the strong coupling constant is evaluated at the scale \i — m/, using the 
NLO RG equation with the appropriate number of active quarks, starting from the input value 
as(M z ) = 0.118. The numerical results for the 0{aa s ) corrections are listed in the first and the 
second columns in Table [2j At the same time, the contribution of the light «-, d-, and s-quark 
to the vacuum polarization is non-perturbative, due to the hadronization effects. In principle, 
it requires a special treatment based on the integration of the experimentally measured spectral 
density within the dispersion relation method (see e.g. |56j), as it was done in [IB]- However, 
this contribution can be estimated by naive use of the perturbative result with effective light 
quark masses. Such estimates are normally in good agreement with the result of the rigorous 
analysis. To estimate the hadronic contribution at two-loop level we use, for the three light 
quarks u, d and s, the value m u = m d = m s = m e ff ~ 180 MeV adopted to describe in the 
lowest order the hadronic contribution to the muon anomalous magnetic moment [ST]. The 
numerical results for the light-quark contribution at KLOE energies are included in the third 
column of Table [2j It is comparable to the contribution of the muon (c.f. Tabled]). Note that 
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Table 2: The s econd- order 0(aa s ) contribution ofb andc quarks, and the light-quark (hadronic) 
contribution to the differential cross section of Bhabha scattering at y/s = 1 GeV in units of 
10 -4 of the Born cross section. 



at 9 = 90° our estimate reproduces the value given in [48J with 20% accuracy, sufficient for 
phenomenological applications at DA$NE. 

Let us now discuss Bhabha scattering at high energies, characteristic to the ILC We consider 
two cases: the Giga-Z option with ^Js = M z \^ and */s = 500 GeV. In the first case, we 
consider the contributions of the leptons and the top quark perturbatively and give an estimate 
of the non-perturbative hadronic contribution due to the five "light" quarks, u, d, s, c and b, 
using the effective quark mass approach. We use the following values for the effective masses: 
m u = irtd = m s = m e ff ~ 66 MeV, m c = 1.25 GeV, and = 4.7 GeV. These values were 
adopted to describe the lowest order contribution to a(M z ) [SB]. The numerical results are 
collected in Table [31 The hadronic contribution is of the same size as the electron vacuum 
polarization contribution and exceeds the one of the muon. Note that our estimate reproduces 
the result of [48] at the scattering angles 90° and 3° with 25% and 10% accuracy, respectively. 
For the small angle scattering, which is of primary interest, our result provides the necessary 
accuracy of 0.1 permille. 

For s/s = 500 GeV, we consider only the contributions of the leptons and the top quark. The 
large electroweak decay width of the top quark serves as an infrared regulator and it suppresses 
the hadronization effects. Therefore, the perturbative result is applicable at energies near 
and above the top-antitop threshold. At the same time, the effective mass approach for the 
contribution of lighter quarks is not reliable at this energy. The corrections become sizable due 
to the logarithmically growing terms and and a rough estimate with an error of about 25%, is 
not accurate enough to match the luminosity precision requirements. Different subsets of the 
second order correction due to the top-quark vacuum polarization to the Bhabha scattering cross 
section are plotted separately in Fig. [81 In this figure we include also the 0(aa s ) corrections 
to the cross section described by the term in Eq. ( 1751) . The numerical results for the two- loop 
corrections due to electron, muon, r-lepton, and top quark vacuum polarizations for a center- 
of-mass energy of = 500 GeV are collected in Table HI Note that the contribution of muon 
and r-lepton can be significantly reduced if one takes into account the corresponding soft pair 
emission, which is justified at this energy from an experimental point of view. 

13 The consideration we make about Giga-Z are valid for LEP1 as well. 
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Figure 8: QED and QCD self-energy ("S"), vertex ('V"), reducible plus one-loop times one- 
loop ("R") and box ("B") contribution to the two-loop top-quark corrections to the differential 
cross section of Bhabha scattering at ^/s = 500 GeV. 
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Table 3: The second-order electron, muon, r-lepton, top-quark, and light-quark contributions 
to the Bhabha scattering differential cross section for the Giga-Z option of the ILC and LEP1 
center- of -mass energy y/s = Mz- The numbers are given in units of 10~ 3 of the Born cross 
section. The top-quark contribution includes the 0(aa s ) term. 
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Table 4: The second-order electron, muon, r-lepton, and top-quark contributions to the differ- 
ential cross section of Bhabha scattering at y/s = 500 GeV in units of 10~ 3 of the Born cross 
section. The top-quark contribution includes the 0(aa s ) term. 



7 Conclusions 

In the present work we derived the two-loop radiative corrections to Bhabha scattering due to 
the vacuum polarization by the virtual heavy-flavor fermion-antifermion pairs. This completes 
the QED analysis of the process at the two-loop level. The result is valid for arbitrary ratio of 
the Mandelstam invariants to the heavy fermion mass, provided all these quantities are large 
compared to the electron mass. The corrections to the cross section are expressed in terms 
of ordinary harmonic polylogarithms and Nielsen polylogarithms of ratios of polynomials in s 
and t invariants. Thus we have complete control over its analytic properties and numerical 
evaluation. The analytical result for the total two-loop QED corrections, which includes also 
the photonic and the electron vacuum polarization contributions [25.36.37J, is now available 

ED]. 

We suggested a new approach which reduces the number of mass scales in the most compli- 
cated part of the calculation. The approach is based on the general properties of the infrared 
and collinear divergencies and could be useful for the high-order perturbative calculations in a 
wide class of processes with a clear mass hierarchy. 

The numerical impact of the perturbative second-order heavy-flavor corrections, including 
the 0(aa s ) contribution for heavy quarks, was studied for the KLOE experiment at DA$NE 
as well as for the GigaZ and the high-energy options of the ILC For the first two applications 
we also provided an estimate of the hadronic vacuum polarization contribution. 

Our result is crucial for the high-precision physics at electron-positron colliders. It removes 
the last piece of pure theoretical uncertainty in luminosity determination at the low-energy 
accelerators and gives the proper account for the top quark effects at the ILC. The accuracy of 
the luminosity determination at low-energy accelerators is now restricted only by the precision 
of the Monte-Carlo event generators for the hard photon and electron-positron pair emission. 
To achieve such an accuracy for the large angle scattering at the ILC one has to perform 
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more careful analysis of the hadronic contribution [48] and take into account also yet unknown 
two-loop electroweak corrections^ 
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A Generalized Harmonic Poly logarithms 

For m e = the two-loop corrections are function of s, t, u, and rrif. It is convenient to use 
the dimensionless variables x and y given by Eq. (Tig]) , and z = R 2 /m 2 = —u/m 2 . To make 
the formulas as compact as possible we introduce six rescaled dimensionless quantities. In 
particular, in the non physical region s < 0, we define: 



P 2 + 4m} -VP 2 s (l-x r f 
x r = / ; ? 2 = ' (^6) 

'P 2 + Am 2 + VP 2 m f X r 



Q 2 + Am 2 -^Q 2 t (1 )a 



2 



Q 2 + Am 2 + m f y r 



R 2 + Am 2 -VR 2 u {l _ Zrf 



,2 



R 2 + Am 2 + VR 2 m f z r 



(77) 



(78) 



'P 2 -4m 2 (l + x h ) 



P 2 + J P 2_ Am 2 m f x b 



IQ 2 - JQ 2 -Am) t a + Vh Y< 

y b = — V =, -L = -^±^L, (so) 



/Q2 + ^Q2 _ 4m 2 m f y b 
vTT 2 -jR 2 -Am 2 u {l + Zb? 



,2 



R 2 + JR 2 -Am 2 m f z b 



14 In the case of e + e~ — » annihilation the two-loop electroweak corrections enhanced by powers of the 

large logarithm ln(M^ vz /s), which are dominant for ^^,500 GeV, were computed in [5S]. This analysis can 
be generalized to the Bhabha scattering by adding the i-channel contribution. 
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The rescaled variable y r is positive for y > 0, while the rescaled variable y b is positive for y > 4 
and is a pure phase for < y < 4. These properties hold also for the variables z, z r , z b and 
The analytical result for the cross section can be expressed in terms of a suitable 
set of GHPLs depending on two of the three variables x, y, and z. Let us consider the pair of 
variables x and y. In this case the set of the weight functions is: 



1 



fo(x) = - 



x 



f-y{ x ) 



\A(4 + x) 
1 



1 



x + y ' " " x + 4 
The GHPLs of weight one are defined as follows: 

G(0; x) = ln(x) , G(a;x) 



/- 



a/x(4 — x) 
1 



x 



(x + y)y/x{A + x) 



dt/ (t) ; 



while GHPLs of higher weight are defined by the iterated integration 



G(a, ■ ■ ■ ; x) 



dtf a {t)G{---;t) 



with the only exception of the weight zero GHPLs, which are defined as follows: 



G(0,---,0;x) 



^-ln n (x). 



(82) 



(83) 



34) 



(85) 



The GHPLs defined in this way satisfy the usual shuffle algebra [IT]. The subset of GHPLs 
which do not involve the weights — y and —y — fi can be expressed in terms of the ordinary 
harmonic polylogarithms (HPLs) of the arguments x r or y r , and weights 1, 0, and —1. The 
table of transformations useful for our calculation reads 



G(0;y 
G(-/j,;y 
G(-4;y 
G(0,0;y 
G(— —fi; y 

G(-4, -fi;y 
G(-fi, -4;y 
G(-4, -A-y 

G(0, 

G(0,0,0;y 



-ff(0;y r )-2ff(l;2/r), (86) 

-H{0;y r ), (87) 

-2 ln(2) + 2#(-l; y r ) - F(0; y r ) , (88) 

H{0, 0; |/ r ) + 2H(0, 1; y P ) + 2#(1, 0; y P ) + 4F(1, 1; y P ) , (89) 

#(0,0;^), (90) 

-C(2) - 2H{-1, 0; y r ) + #(0, 0; j/ r ) , (91) 

C(2) + 2 ln(2)#(0; y r ) - 2//(0, -1; y r ) + #(0, 0; y r ) , (92) 
2 ln 2 (2) - 4 ln(2)//(-l; y r ) + 4//(-l, -1; y r ) - 2H{-1, 0; y r ) 

+2 ln(2)F(0; y r ) - 2/7(0, -1; y r ) + (0, 0; y r ) , (93) 

2C(2) + jy(0, 0; j/ r ) + 2^(1, 0; j/ r ) , (94) 
-//(0, 0, 0; y r ) - 2//(0, 0, 1; y r ) - 2/7(0, 1, 0; y r ) - 4if (0, 1, 1; y r ) 
-2/7(1, 0, 0; y r ) - 4//(l, 0, 1; y r ) - 4//(l, 1, 0; y r ) 

-8//(l,l,l;y r ), (95) 

-H(0,0,0;y r ), (96) 
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G(-fi, -4, -//; y 
G(-4, -4, -//; y 

G(0, 0,-^2/ 

G?(-/i, 0, -//; y 
G(0,-n,-n;y 

C(-4, -/i, 

G(0;x 

G(0,0;x 
G(/i, //; x 
0; x 

G(0,0,0;x 



/i, 0; x) 



3C(3) + C(2)i?(0; y r ) + 2/7(0, -1, 0; y r ) - H(0, 0, 0; y r ) , (97) 
2C(3) - 2C(2)fT(-l; y r ) - 4#(-l, -1, 0; y r ) + 2#(-l, 0, 0; y r ) 

+((2)H(0; y r ) + 2H(0, -1, 0; y r ) - H(0, 0, 0; y r ) , (98) 
2C(3) - 2C(2)iJ(0; y r ) - H(0, 0, 0; y r ) - 2H(0, 1, 0; y r ) 

-4C(2)F(1; y r ) - 2F(1, 0, 0; y r ) - 4H(1, 1, 0; j/ r ) , (99) 

-4C(3) - 2C(2)#(0; y r ) - #(0, 0, 0; y r ) - 2#(0, 1, 0; y r ) , (100) 

2C(3) - H(0, 0, 0; y r ) - 2H{1, 0, 0; y P ) , (101) 

-§C(3) + 2H(-1, 0, 0; Vr ) - if (0, 0, 0; y r ) , (102) 

2ff(-l;x 6 ) -if(0;x 6 ), (103) 

vr + zif (0; x b ) , (104) 

4/f (-1, -1; x b ) - 2H(-1, 0; x b ) - 2/7(0, -1; x 6 ) + ff(0, 0; x 6 ) , (105) 

-3C(2) -ff(0,0;x 6 ), (106) 

-iC(2) + 2zif (0, -1; x 6 ) - zif (0, 0; x 6 ) , (107) 
8/f (-1, -1,-1; x b ) - 4if (-1, -1, 0; x b ) - 4ff (-1, 0, -1; x b ) 
+2H(-1, 0, 0; x b ) - 4ff(0, -1, -1; x b ) + 2/f (0, -1, 0; x b ) 

+2H(0, 0, -1; x 6 ) - if (0, 0, 0; x b ) , (108) 

-2C(3) + C(2)#(0; x 6 ) - 2ff (0, 0, -1; x b ) + if (0, 0, 0; x 6 ) . (109) 



The auxiliary two- loop box 5-functions listed in Appendix iBlinvolve three GHPLs which depend 
on two different kinematical variables. These GHPLs can be expressed in terms of the logarithms 
and Nielsen's polylogarithms depending on the rescaled variables: 



G(-x - fi;y) 



G(-x, -n;y) 



dw- 



1 



(w + x)yjw(4 + w) ' 

(ln(x 6 + y r ) - \a{x b y r + 1)) 



(1 - x b )(x b + 1) 

dw 

n w + xj 



(110) 



dr 



V r ( 4 + 



G(-x,Q,-fi;y) 



\n 2 (y r ) - ln(x 6 + y r ) \n(y r ) - \n(x b y r + 1) hi(y r ) + Li 2 

V. Vr 

-Li 2 {-x b y r ) , 

y i r w i r q i 

dw / dq— / dr — -^=^= , 

o w + xj qj ^r(4 + r) 

ln(l + x 6 )j -21n(y r )[ln(x 6 y r + 1) + ln(x 6 + j/ r ) - \n(y r )} - 8((2) 



(111) 
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x b y r + l 
x b y r 



-2Li 9 



+2Li 2 [ XbV ; + l J + 2Li 2 



Xfc + y r 

x b +y r 



■2U 2 [ X -^±l\-2uJ Xb+V 



1 + x b J \1 + x b 

, , . . . — w . , 21n(x 6 )m(y r )+41n(y r )ln(l-y, 

y r (l+X 6 )/ VT/r-Cl+Xfo)/ 

+4Li 2 (l -y r ) + 8Li 2 (y r ) +4 ln 2 (l +x 6 ) ln(y r )l + [ln(x 6 +y r ) -ln(y r )][2C(2) 
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+2Li 2 (^±^U2Li 2 (^±± 

x b y r + l 



2Li 2 



Xb + Vr 

y r (l + x b ) 



+ ln(a; 6 |/ r + l)2C(2) + 2Li 



x b y r 



+2 ln(x b ) ln(y r ) + - In 2 (y r ) - 2Li 2 (y r ] 



-2Li, 



l + x b 
x b y r + l 



-2Li 3 P*±J 
V l+x b 



-2Li3|^ l+ 2L i3 



2W, 

y r (l+a; 6 ) 



-2Li, 



2Li 2 (f/ T 
y r (l+x 6 ) 



+2Li 3 



^(l + Xft) 



-In(j/ r 
-U 2 (-x b y 



2Li2 (3*+^ _ 2Li2 Z^+jA _ 2Lij (jMtlL) _ Li2 (_» 

x b y r J \l+x b J \y r (l+x b )J V x b 



+ l\n 2 (y r )[3\n(x b +y r )-4\n(y r )-6\n(x b )] + 2U 3 ' 1 
o 



+2Li 3 (l + x fe )-Li 3 ( J+Li 3 f -— )-Li 3 (-x 6 y r )+Li 3 (-x 6 ), (112) 



where the Nielsen's polylogarithms are related to HPLs as follows: 

Li n (a) = H(0, • • • , 0, 1; a) , S n>m (a) = if(0, • • • , 0, 1, • • • , 1; a) . 



(113) 



A.l Analytical Continuation 

The result for the two-loop corrections in Section [5] is expressed in terms of the auxiliary 
functions which are given in Appendix [B] in the non-physical region s < 0. The corresponding 
expressions in the physical region s > can be obtained by analytical continuation to the 
complex value of P 2 : 

P 2 = -s- ie, e -> 0+ . (114) 

Not that in the physical region Q 2 and R 2 are real and positive. Let us consider the analytical 
structure of the rescaled variables, Eqs. fl76H81~]) . The variables y r and z r are positive and 
vary from to 1 for < Q 2 < oo and < R 2 < oo. The variable x r varies from to 1 for 
< P 2 < oo. For positive s the variable x r becomes complex. In the region < s < 4m 2 it is 
a pure phase: 



'4m 2 f — s + i^/s 

X r = V , = (P* , (115) 

4m? — s — iy/s 



where 

In the region s > 4m 2 we have 
with 



arctan , / = . (116) 

1 ' Amj — s 



x r = —x' r + ie , (H7) 



y/s — 4 s — Am 2 r 



f 



:n8) 



y/s + * Is — Am 
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The HPLs of the variables y r and always real, while the HPLs of Xf £1X6 complex for s > 0. 

In particular, their imaginary part in the region above the heavy-flavor threshold, s > Am 2 , is 
defined when the analytical continuation of the logarithm is specified: 

H(0- x r ) -> H(0- -x' r + ie) = H(0; x' r ) + in . (119) 

The case of the variables x b , y b , and z b is more complicated. The variables y b and z b can be 
complex. For < Q 2 , R 2 < Am 2 , y b and z b are pure phases. Note that the expressions of t and 
u in terms of y b and z b are invariant under the inversion 

Vb -> - , ^ -> - • (120) 

2/6 ^6 

This means that for < Q 2 ,R 2 < 4m 2 it does not matter whether we give to Q 2 and R 2 a 
positive or a negative imaginary part. Let us choose Q 2 —> Q 2 — ie and R 2 — > R 2 — ie. Then 
we have 



V/Q2 + i Urn? - Q 2 

y b =— V , =e a +, (121) 

y^Q 2 - iJAm) - Q 2 



R 2 + ijAm 2 - R 2 



z b = \ ' = e 12 ^ , (122) 



R? - ijim 2 - R 2 



with 



Um 2 t 

$ = facte* (123) 



Um 2 

£ = arctan W - 1 . (124) 

When Q 2 and R 2 vary from 4mj to infinity, y b and Zft are real and positive and vary from 1 to 
0. The variable x b is a pure phase for < P 2 < Am 2 . It varies from 1 and when P 2 varies 
from Am 2 to +oo. In the physical region we have 

x b — > —x' h + ie, (125) 

where 



1 s + 4m 2 f — y/s 

x' b = v , , (126) 

's + 4mj + y/s 



i.e. x b is negative and varies from -1 and when s varies from to oo. The HPLs of x b are 
complex in the physical region. Their imaginary part is defined in the same way as in Eq. fll 19j) . 
Finally, the analytical continuation of the three GHPLs of Eqs. ( Ill OH 112"]) are defined by the 
analytical properties of the functions Li 2 and Li 3 
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A. 2 Mellin-Barnes Expansion of the GHPLs 

To study the low and high energy behavior of the two-loop corrections one needs the expansion 
of three GHPLs of two kinematical invariants, Eqs. (I11CW112]) . in the limits s ^> m'j and s -C m 2 . 
To perform the expansions we apply the inverse Mellin-Barnes transformation to the integral 
representation of the GHPLs. In this Section we describe the technique and present the result 
of the expansion. We start with the GHPL of weight one in Eq. (11101) : 

G{—x — fj,; y) = f dw- ) (127) 

Jo [w + x)^/w{4 + w) 

By changing the integration variable w = yr we obtain the following integral representation: 

If 1 !/ y x-i/ 4^- 1/2 



G(-x-ii;y) = - dr- 1 + -r 1 + — • (128) 
x J r V x J \ yr J 

Then we apply the inverse Mellin-Barnes transformation to the square root in the integrand: 

4 V 1/2 1 1 r°° (\ \ ( 4 V 

' daT[=- + a)T(-a) — . (129) 



- 



yrj 2mT(l/2)J_ loo \2 J v ' \yr 

The integrand has two infinite series of poles in the a complex plane: 

"left hand side poles" at a = —n — - , for n — 0, 1, 2, • - • ; 
"right hand side poles" at a = n , for n = 0, 1, 2, ■ ■ • . (130) 

One has to choose the integration contour in Eq. f)129p so that 

-i<Re((x)<0. (131) 

Then we find 

G <-* - " : v) = hfm \ II iaT S + °) r (3 ' ' (132) 

The integral over r can be evaluated analytically 



I(a)= drr-^h+^r) = 2 Fi 1, -a, 1 - a; ), (133) 

Jq \ X / (J \ X/ 

where 2F1 is the hypergeometric function. In order to obtain the asymptotic expansion in the 
y — ► 00 limit at fixed r = y/x it is sufficient to close the integration contour in Eq. (11321) on 
the r.h.s. of the a complex plane. One finds 

00 

G(-x-^y) = J29 { n' L) (r,y)(m n+1 > (134) 

n=0 
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where the functions gn'^ can be easily obtained e.g. with the help of [13]. For n — 0, 1, 2 
they read 

g^ L \r,y) = -r(ln(r + l)-ln(j/)) , (135) 

£i 1,L) (r,2/) = -2r[(ln(r + l)-ln(y))r + r-l] , (136) 

9^\r,y) = -r [6(log(r + l)-log(y))r 2 + (7r-6)r + 3] . (137) 

To obtain the asymptotic expansion of this GHPL in the limit y — > for fixed r it is necessary 
to close the integration contour in Eq. (11321) on the left hand side (l.h.s.) of the o complex 
plane. One obtains 

oo 

G(-x - fi; y) = Y^9 { n' S) (r,y)y n , (138) 

n=0 

(1 S) 

where the functions gn ' for n = 0, 1, 2 read 



(i.s) 
0o 



/- arctan (y/r) , (139) 
y 

g™ = "^^^ " ^ , (140) 



(i,s) y/r(r-3) + 3 arctan (y^) 41 
~ 128 

The expansion of the remaining two GHPLs is almost identical. By inverse Mellin-Barnes 
transformation the GHPL of weight two in Eq. (II 111) can be written as 

rv \ r w 1 
G(—x,—(j,;y) = / dw / dr- 



w + x J ' ^r(4 + r) 

y i /-i i / 4 

/ <is- 1 H 

w + r Jq s \ ws 

1 1 [ iM , T(Ua) T{-cr) n fv 1 

' da^ J — 4 / ciw w" 



27rzr(l/2) J_ ioo a J w + x 

i i i r°° r (± + a) r(-o-) 



2nr(l/2)xj_ 400 ct(1-<t) 
x 2 F 1 (l,l-<7,2-<7;-£) . (142) 

Again, the integration contour has to satisfy Eq. (I13ip . By closing the integration contour on 
the r.h.s. of the complex a plane, in the limit y — > oo one finds 

oo 

G{-x-^y) = X>?' L) (T,y)(l/</) n . (143) 

n=0 

The functions gn'^ for n = 0, 1, 2 read 

^ 2 ' L) (r,y) = ln(y)ln(r + l) + Li 2 (-r), (144) 
g^ L \r,y) = 2r[ln(y) + ln(r + 1)] , (145) 
g?> L \r,y) = r[3(\n(y) - ln(r + l))r - 2r + 3] . (146) 
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Alternatively, by closing the integration contour on the l.h.s. of the complex a plane, which 
corresponds to the limit y — > 0, one obtains 

oo 

G(-x,-Ai;y) = $>i 2 ' S) (r, y)y n , (147) 



ra=0 



with 



^ 2 ' S) = 2y|(v^-arctan(v^)) , (148) 

(2,s) V? ( V^( r - 3) + 3 arctan (yg)) 

^ " 36r ' li49j 



(2,s) = Vf (y^(r(3r - 5) + 15) - 15 arctan (JF)) 
92 ~~ 1600r 2 ' 1 J 



For the GHPL of weight three in Eq. fillip we have 

ry 1 T w 1 
G(—x, 0, —fx; y) = I dw / dq- I dr- 



w + x J q J y/r(4 + r) 

y 1 f w 1 f 1 1 f 4 
dw / ds- / dq- 1 H 

o w + x Jo s Jo q V S( l 



doT ( a + - ) r(- ( x)4 CT / ^ du>^^- x 



2m r(l/2) 7_ ioo V 2) y J Jo w + x 



i i r r (a + ±) r(-<r) 



27rzr(l/2)x7_ ioo a 2 (l-a) 
x 2 Fi(l,l-a,2-cr;-|) . (151) 

Once again, the condition in Eq. (I13ip has to be satisfied by the integration contour in Eq. (I15ip . 
By closing the integration contour on the r. h.s. of the complex a plane one obtains the 
asymptotic expansion in the limit y — > oo: 

oo 

G(-x,0,-Ky) = J29 { n' L) (r,y)(l/y) n . (152) 

n=0 

The first three gn' L ^ functions are given by 

9^ L \r,y) = i(31n 2 (y) + 27r 2 )ln(r + l) + ln(y)Li 2 (-r)-Li 3 (-r), (153) 
g[ 3 ' L \r,y) = -2r[ln(y) - ln(r + 1) + 1] , (154) 
9i 3,L \r,y) = ~r[6(ln(j/)-lB(r + l))r-r + 6]. (155) 
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By closing the integration contour in Eq. (11511) on the l.h.s. of the complex o plane one finds 
the asymptotic behavior of the GHPL in the y — > limit: 



G(-x,0,-KV) = J2g { n' S) (r,y)y n , (156) 

71=0 

where the functions gn' for n — 0, 1, 2 read 

9 ( o' S) = 4y|(v^-arctan(v^)) , (157) 

(3,s) \/^(Vt(t-3) + 3arctan( v /7)) 

* = 54^ ' ( 158 ) 

(3,s) (y^(r(3r - 5) + 15) - 15 arctan (yg)) 

^ = 4000^ • (159) 

It is interesting that, though the expansion of the three GHPLs discussed above in the s -C mj 
limit involves arctan(y / r) terms, these terms completely cancel in the expansion of the auxiliary 
functions Bi, B 2 , and B 3 . 



B Auxiliary Functions 

In this Appendix we collect the expressions for the auxiliary functions used in the paper which 
are valid in the non-physical region s < (P 2 > 0). The analytical continuation to s > is 
discussed in Section IA.1I The dimensionless variables x and y, used in the explicit formulas 
below, are related to s and t via Eq. (fl9|) . 

B.l One-Loop Functions 

Vacuum polarization: 

i 

TlS l \s) = Y,n<i Wl \s)(D-4y + 0{(D-4f) , (160) 



i=0 



where: 



tt(".°) 5 4 1 x 2 + 2x-8 

9 3x 3 X y/x(x + 4) 



n ( 1/ . 1 ) 



Dirac form factor: 



14 16 1 



(4x-64 + 5x 2 )G(- / u;x) 



27 9a; I8xy/x{x + 4) 

-3(2x - 8 + x 2 )G(-4, -/i; x) . (162) 



o 

^(s) = Fl ll>i) (s)(D - Af + 0(D - 4) , (163) 



i=-l 
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where: 

*f - 1 - ln (J) • < 164 > 

ff •»> = -1 + i C (2) - I G (0; *) + f 5 + i G (0; *)) In (*) - ? In 2 f . (165) 



e 



Box £>-functions: 

o 

( s ' = £ 5 f ' l) (M) - 4)< + 0(I> - 4) , J = 1,2,3, (166) 



i=-2 

where: 

B (u,-2) = 16(x + y) 2 
1/ 

^ = 8( " 2 + 3:?/ + y2) H- 8(3: + y)2 G(0;x), (168) 

y y 

3(11,0) = 2(2sy + 8a; 2 C(2) + 10^(2) + 5y 2 C(2)) 

1 1/ 

+ 2(2x +xy + y ) ^ + 2 + ^ + 2 ^ + Q _ ^ 

1/ 

+ 2(2x2 + ^ + 1/2) (G(0; x)G(0; y) - G(0, 0; y)) , (169) 

(170) 

= 8y+ 8{2x2 + 2xy + y2) G(0; X ), (171) 
= _ 2(16^ + loly + 5^) C (2) + 2(4^ + 2,y + y>) 

y y 

-2{x - y)G(0; x) + 2{x + y)G(0; y) + 2(2a; + y)G(0, 0; x) , (172) 
Bl ll >-V = -1^!, (173) 

Jf" 1 ) = - 8( " 2 + ^ + ^-^(0;x), (174) 

bT> - - 4(aCy + ya - 4gaC(2)) -^(G(0; g ) g (0;y)- g (0,0;y)) 
_ 4(^ + ,y + y 2 ) 

y 

B.2 Two-Loop Functions 

Vacuum polarization: 

nf ) ( s ) = n^' 0) ( s ) + o( J D-4), (176) 
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where: 



n 



(2i,0) 



x + 4 



3(-6 + x)G(-At; x) + 4(2 - x)(2G(-4, -//; x) + G(0, -fi; x))] 



12xa/x(x + 4) 

+ 24^2 52x ~~ 5x2 ~~ 8 ( 7 ~ 2x ~ 3a;2 ) G (-/ i ? -A 4 ; 

+ (4 - x 2 )(16G(-4, -n, -fx; x) + 8G(0, -li, x) - 16G(-/i, -4, -//; x) 



-8G(-li,0,-li;x)) 
Dirac form factor: 



where: 



(21,0) 



Ff a \s) = F? a ' \s) + O(D-4), 
: [(x-4)(46- 19x)G(//,0;x)] - 



(177) 
(178) 



36x^(4 - x) LV ~' v " " ' /j 1296x2 

-x(3960 - 1590x)G(0; x) + (1296 - 216x 2 )G(/z, //, 0; x)] , 



[8568x - 3355x 2 



Box £>-functions: 



(179) 



(180) 



i=-2 



where: 



5| 



(21,-2) 



8(x + y) 2 (5y-12) , 8 (y-2)( 2/ + 4)(x + y) 2 



G(-n;y) , 



(21,-1) 



V 3yV^ + 4 ) 

4(204x 2 + 444xy - 13x 2 y + 204y 2 - 41xy 2 - 13y 3 ) 

2V 

8^ + y) 2 , 4(x + y) 2 (5y- 12) 

-yu; y) — G(0; x) 



(181) 



3y 
1 



V 

4(y + 4)(-34x 2 - 74xy + 2x 2 y - 34y 2 + 7xy 2 + 2y 3 



4(|/-2)(y + 4)(x + y) 2 

3y 2 



V 



(3G(-4, - /i;2/ ) + G(0;x)G(- /U ; 2 /)) , 



(182) 



^1 



(2J.0) 



2(x - 4)x(y - 2)(y + 4)(x + y) ; 
3y 2 v / y(y + 4)y/x(4-x) 



G(li, 0; x)G(-x - ii\ y) 



Vv(y + 4 ) 



2(y-2)(y + 4)(a; + y) ; 
3y 2 



G(-x,0,-ii;y) 



4(y - 2)(y + 4)(x + y) 2 y + 4 2 

-G(-/i, -//, -/i;y) + T ^ 7 ^(-1736x - 4732xy 



135y 2 



+82x 2 y - 1976y 2 - 16xy 2 + 9x 2 y 2 - 158y 3 + 18xy 3 + 9y 4 + 180x 2 C(2) 
+360xyC(2) - 90x 2 y((2) + 180y 2 C(2) - 180xy 2 C(2) - 9(h/ 3 C(2))G(-/i; y) 

- 2(y Q + 2 4) (-34x 2 - 74xy + 2x 2 y - 34y 2 + 7xy 2 + 2y 3 )G(-4, -//; y) 
3y 2 
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2(y-2){y + 4)(x + yy 
3y 2 



(3G(-4, -//; y)G(0; x) + 9G(-4, -4, -//; y) 



-G(0; -//; y)) - ^^(4x 2 +4xy-5x 2 y-4y 2 -5xy 2 -4y 3 )G(0, -fi; y) 

2<y4 + 2 V \ -A2x 2 - 82xy + 12x 2 y - 26y 2 + 17xy 2 + 10y 3 )G(0; x)G{-fi; y) 



+ 



9y 2 

2(y-2)(y+4)(x + y) 2 {G{Q _ x)G{Q ^. y) g,-//; y)+G(0, 0; y))] 

2(x - 4)x(12x + 20y - 5xy + 3y 2 ) 1 AA 9 9 

+ - „ , . = — G A*, 0; x - ^(144x 2 + 288xy - \\x 2 y 

9y 2 A/a;(4 - x) 9y^ 

+144y 2 - 40*y 2 - 17y 3 )G(-/i, - m y) - ^E±JH g(-^ -4, - m y) 

y 

2(2x 2 + 2xy + y 2 ) 2{2x + y) (n(n 

G(-n, 0, -n\ y) (G(0; x)G(-n, -//; y) 



3y 

+G'(„. /i, 0; *)) + 2(2X2 + ^ + %2) G(0, //. y) + J~ 2 (-180x 2 - 372xy 



3y 



27y 2 



+41:r 2 y - 156y 2 + 12Lry 2 + 65y 3 )G(0; x) - ^^^G(0, 0; x) 

yy 

(x + y) ^ + _|_^ 4488;r 2 + 120 36a;y - 1082x 2 y + 4488y 2 - 2359xy 2 



15 ~ ' 405y 2 

-1082y 3 - 540x 2 C(2) - 1080a;yC(2) + 225:r 2 yC(2) - 540y 2 C(2) + 450xy 2 C(2) 
+225y 3 C(2)), (183) 



8(2, 2 + 2,y + y 2 ) ^ _ % + 3fa,-2)( y + 4) 



9y 2 



(184) 



4(480x 2 + 480a;y - 56x 2 y + 204y 2 - 56xy 2 - 13y 3 ) 4(5y - 12) , 9 
(2x + 2a;y 



+y 2 )G(0;x) 



27y 2 

8(2a; 2 + 2xy + y 2 ) 
% 



9y 2 



G(-//, -/i;y) 



y + 4 



Vv(y + 4 ) 



— -(40x 2 + 40:ry 
9y 2 



r 2 ±17 2 r; 2 3^/ \ , % - 2) (2a; 2 + 2xy + y 2 ) 
-5xy+17y - 5xy - y d )G(-/i; y) + — (3G(-4, -//; y) 



+G(0;x)G(-v,y)) 
2(x - A)x(y - 2)(y + 4) (2a; 2 + 2xy + y 2 ) 



3y 2 v / y(y + 4)^(4 -rr) 



3y 2 



0; x)G(-x - n\ y) 



(185) 



y + 4 



2(y-2)(2x 2 + 2a;y + y 2 



0, -n; y) - 6G(-fi, -/i, -fi; y)) 



Vv(y + 4 ) 

^- (5872a; 2 + 6112xy - 464x 2 y + 1976y 2 - 224xy 2 - 18a; 2 y 2 + 158y 3 

135y^ 

-18xy 3 - 9y 4 - 360x 2 C(2) - 360a;yC(2) + 180x 2 yC(2) - 180y 2 C(2) 

4 
3y 



+ 180xy 2 C(2) + 90y 3 C(2))G(- / u; y) + 7^(40x 2 + 40xy - 5x 2 y + 17y 2 - 5xy 2 
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3y 2 



+9G(-4, -4, -n- y) - G(0; x)G(-x, -//; y) + G(0; x)G(0, -//; y) 

-G(0, 0, -/x; y) + G(0, 0; y)) - -^(8x 2 + 4xy-10x 2 y-4y 2 -5xy 2 

9y l 

3 4(48x 2 +44xy-15a; 2 y+13y 2 -10a:y 2 -5y 3 ) i 
-4y J )G(0, -fi; y) + — G(0; x)G{-^ y) 

2(x - 4)x(24x + 16y - 10a;y + 3y 2 ) ^. 1 . 9 _ 9 

+ — ; = y —GU, 0; x) -(288a; 2 + 288xy - 46a; 2 y 

9y 2 



9y 2 \Jx(4 — x) 

+U4y 2 - 52xy 2 - 17y 3 )G(-/i, -fi; y) - (2a; 2 + 2xy + y 2 ) (-G(-/i, -4, -//; y) 



+ 1 G(0 ; B ) + %^G(0, 0; ,)) - !M±|W) G( _,, o, _„ y) 



15 



9y 2 



3y 



-^(2x + y)(G(/i, fi, 0; x) + G(0; x)G(-fi, -fi; y)) - ^p(^ 2 + 432a;y 

-112a; 2 y+ 156y 2 - 136a;y 2 -65y 3 )G(0; x) + ^!±^±^!i G(0 ) 

3y 

fTT^r(16176a; 2 + 16896xy-3004x 2 y+4488y 2 -3004a:y 2 -1082y 3 -1080a; 2 C(2) 

(186) 



405y 2 



-1080xyC(2) + 450x 2 yC(2) - 540y 2 C(2) + 450a;y 2 C(2) + 225y 3 C(2)) , 
8a; 2 (5y-12) 8a; 2 (y-2)(y + 4) 

V 3yVl/(i/ + 4) 

4(204a; 2 - 36xy - 13a; 2 y - 36y 2 + 15a;y 2 + 15y 3 ) 4a; 2 (5y - 12) 



(187) 



27y 2 
y + 4 



+ 



9y 2 



G(0;x) 



9y 



(34a; 2 — 6xy — 2x 2 y — 6y 2 + 3xy 2 



+-^-G{-fi,-fi;y)- 

3 y Vv(y + 4 ) 

+3y 3 )G(-/i; y) + J (3G(-4, y) + G(0; x)G(-|x; y)) 

2(a;-4)a; 3 (y-2)(y + 4) 

3y 2 y/y(y + 4)^(4 -x) 



(188) 



G(/x, 0; x)G(-x - n\ y) 



y + 4 



2a; 2 (y — 2) 

— 0, -/i; y) - 6G(-/i, -//, -//; y)) 



Vy(y + 4) 

^r( 1736a; 2 - 1020a;y - 82a; 2 y - 1020y 2 + 60xy 2 - 9a; 2 y 2 + 60y 3 - 180a; 2 C(2) 

\3hy z 

+90x 2 y((2))G(-i2; y) + -^(34a; 2 - Qxy - 2x 2 y - Qy 2 + 3xy 2 + 3y 3 )G(-4, -fi; y) 



2x 2 (y-2) 
3y 2 



3y 

(G(0; x)G(-x, -k y) - G(0; x)G(0, y) - 3G(0; x)G(-4, -/x; y) 



-9G(-4, -4, -//; y) + G(0, 0, -//; y) - G(0, 0; x)G(-/x; y)) 



4a; 2 (5y-4) 
V 



G(0,-/i;y)4 



2 ( 1 4a; 2 - 2xy - 4a; 2 y - 2y 2 + xy 2 + y 3 ) 
3? 



G(0;x)G(-v,y) 
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2(x — 4)x(12x — Ay — 5xy) 
9y 2 ^x(A-x) 



G(p, 0; x) 



405y 



2 (4488x 2 - 23A0xy - 1082x 2 y 



1 



-3060y 2 + 195xy 2 + 195y 3 - 540x 2 C(2) + 225x 2 yC(2)) + -^(144x 2 - llx 2 y 

V 

4x 2 

-12xy 2 + 12y 3 )G(-p, -p; y) + — (3G(-p, -4, -/i; y) + G(-/x, 0, -p; y) 

2(180x 2 -12xy-41x 2 y-36y 2 + 15xy 2 + 15y 3 ) ^ n 
-G(0, -/i, -/i; y)) + — -= G(0; x) 



27y 2 



(189) 



C Expansions of the Cross Section 

The non-logarithmic part of the second order correction in Eq. (Q can be written as follows 

,2\ / r \2 



S, 



(2) 



-2 



1 + ln 



In ( + QjiV c ( T -| T ^J /(p,0 



(190) 



where the first term is determined by the soft emission and /(p, £) is a function of two dimen- 
sionless variables: p = m 2 /s and £ = —t/s. The small-mass expansion of the function /(p, £) 
(p = m 2 /s, £ = —t/s) is of the following form 



(191) 



n=0 



where / n (p, £) depend on p only logarithmically. For the leading term we obtain 



(e 2 -e+i) 2 (i 

£ 2 19 



/o(p, = sr^" { ~ In 3 (P) + In 2 (p) 



27 +3 WO 



+ ^ (p) 

32£ 4 - 46£ 3 + 33£ 2 + 8£ - 4 

+ W 
20£ 4 -31£ 3 + 60£ 2 -31£ + 20 



40 T . / j-\ H651 , , , 

+ Y Ll2(0 + ^r/- ln(p) 



C(2) _ - g + 1) y - 7 g + 4) 



18£ 2 

8£ 4 - £ 3 - 15£ 2 + 17£ - 4 



ln(l-0 + 



6£ 2 

20£ 4 - 67£ 3 + 141£ 2 - 1 1 2£ + 74 



18£ 2 



MO 



12£ 2 

(2£-l)(£ 2 -£+l) 



ln(0 2 - 



(2£-l)(4£ 3 -3£ 2 + 4) 



6£ 2 



ln(01n(l-0 



3£ 



C(3) 



- I) 2 (e 2 - e + 1) 
9£ 2 



ln 3 (l-0 



196£ 4 -311£ 3 + 258£ 2 + 13£-38 c(2) _ 2 (2£ 4 -9£ 3 + 16£ 2 -ll£+4) ^ ^ (2) 



18£ 2 

12£ 4 -20£ 3 -£ 2 + 24£-4 



3£ 5 



6£ 2 



ln(0C(2) + 2(1 eH 3 f 2 e + ln(l - 0Li 2 (0 
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7 (16e 4 -23e 3 +48^ 2 -23^+16) 1 M 20£ 4 -58£ 3 + 81£ 2 -58£ + 20 , 2/1 
H — - — — ztttt^ -ln(l-£)H ln^(l— ^) 



54£ 2 



18£ 2 



4^-6^ + 7^-4^^^ 2 40£ 4 -50O + 9O + 37«£-20 



ln(0ln 2 (l-£) + 



12^ 

e 4 -3e 3 + 4e 2 -e + i 
3^ 2 

20£ 4 + 8£ 3 - 84£ 2 + 92£ - 55 
W 2 



18£ 2 



ln 2 (0ln(l-0 + 



4^4 _ 2 ^3 _ 22^2 + 3i^_4 
36^ 



HO H^-0 

ln 3 (0 



ln 2 (e) 



(^ 2 -^ + l)(2e 2 - 7^ + 12) 



3? 



HO^HO 



112£ 4 - 449£ 3 + 1011£ 2 - 836£ + 562 

5ie 



HO 



2d - a (e - e+i) Li3(1 _ + (^+i)(^+4) Li3(0 



W - 1) 



3£ 2 

e 2 



3e 

(i-e+e 2 )(^-2C(3)+^in( P ) 



2-£ 



WO 



(192) 



in agreement with the result of Refs. [27|4"5] . We observe that the functions /$ depend on the 
charge of heavy fermion through the contribution of the two- loop irreducible self-energy 
diagrams, which are proportional to Q 4 , while all the other graphs that we consider in the 

2 rpi 



present work are proportional to Q 2 f . The next-to-leading term is new and reads 



A(p,0 



m-D ( ^+d (2^+2) [ln2 (p)+4Li2(0+12] + ^{ (e _ 1)(2 ^3 



+2£ 2 ) ln(l - £)] 



2(12£ 4 



+5£ 2 - 5£ + 2) + 2(£ 2 - £ + 1)[(4 - 2£ + e - 2£ 3 ) ln(£) + (£ - 1)(2 + £ 

4r(40£ 5 -54£ 4 +50£ 3 -i7£ 2 -i2£+8) + 

^ 2^ 
-5£ 3 + 13£ 2 - 5£ + 12) + (8 - 6£ + 9£ 2 - 6£ 3 + 8£ 4 ) ln(l - £)] H 1 - 
1 



20 



2(12£ & - 21£ 4 + 26£ 3 - 26£ 2 + 21£ - 14) - 2(4 - 6£ - £ 2 + 8£ 3 - lOf 



+80) ln(l - + (8 - 15C + 12£ 2 - - 70 + 8£ 5 ) HO HO 



(Q} ~ i) fs 



(2-3£+4£ 2 -4£ 3 + 3£ 4 -2£ 5 )ln(p) 

_( 2 _ 3 £ + 3£ 2 -£ 3 ) ln(£) 
The expansion in the large-mass limit takes the form 

oo 



(193) 



(194) 



n=0 



where the leading n = term vanishes because of the renormalization condition and f n (p,x) 
depend on p only logarithmically. For the next-to-leading term we obtain 



Mp,0 



955£ 3 - 3926£ 2 + 3926£ - 955 12£ 3 - 19£ 2 + 14£ - 6 



1350£ 



C(2) 



3e 3 +e 2 -e-3 



In (1-0 
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af-Bf + Be-2 , 5f-22f + 22C-5 

20£ V ^ J 30£ ^ 

20f 3 -78f 2 + 93f-58, _ , M 12f 3 - 19£ 2 + 14f - 6 . . 

■ In (0 + In (1-0— " " ln £ 



90£ v>/ v w 30£ 

+l He + z + 4) in 2 (o + 4(e3 ~ 2 ^ +2e " 1} Li 2 (o 

Finally, the next-to-next-to-leading order in the s <C mj expansion reads 
177763^ - 405359O + 6761940 - 405359£ + 177763 



/ 2 (P,0 



2116800£ 

3 (40 - 170 + 160 - 12f + 2) 20 + 150 + 60 + 15£ + 2 



C(2) + -2—! s s ^— ln(l - £) 

280£ SV ; 840£ V ^ 

(0-4^+1) (20 -3^ + 2) _ 53Q-14lO + 222Q-141^+53 



560£ v 1680£ 

330-1390 + 2130-212^ + 52 _ 4Q-17Q + 16Q-12e + 2 _ 

1680£ nl ^ j 280£ nl ^ j 111 U 

+ 5^0 ^ " ^ ~ 2 ^ + l) ln( ° 2 ~ (e2 ~ 4e + 3 5 f ~ e + 1) Li2 ^ 

+ (^-l) 449(1 - 5 ^ + 6 ^- 5 ^ + e4) . (196) 
w y 10800£ V ; 
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